STELLAB and OMEGA (Stellar yields + Faint Supernovae)

Documented by Jacob Brazier.

Note This notebooks require an experimental yields table, which is not part of NuPyCEE.

See Côté et al. (2016) and the OMEGA/STELLAB/SYGMA notebooks for further information.

Faint supernovae are low luminosity core-collapse supernovae characterised by their relatively small generation of Nickel-56 Nomoto et al. 2013. Nuclei that originate from the silicon burning regions of a star, undergo fallback if the explosion energy is too low to overcome gravitational attraction. The mass of the content that undergoes fallback is called the 'mass cut'. Anything above this mass cut is pushed outwards by the explosion, while the material underneath undergoes neutronisation/destruction by the formation of a black hole. The material that is ejected tend to include oxygen and carbon nuclei. Mixing may also occur at the masscut, so that near-iron nuclei may still be ejected. The mixing and fallback model is described in more detail by Umeda & Nomoto 2002

This project uses OMEGA (One-zone Model for the Evolution of Galaxies) to simulate the chemical evolution of galaxies and to plot the expected abundance ratios of titanium, oxygen, iron, hydrogen and carbon. The main plots show the how the abundance of these elements (with respect to hydrogen) change with both time and stellar metallicity. Fe/H is assumed to be a measure of metallicity, as higher Fe/H values indicate a higher proportion of metal relative to hydrogen. Fe/H is also a rough indicator of the age of stellar populations as metallicity tends to increase with time, and the age of stars can be difficult to determine. We plot the predictions of OMEGA with and without the yields of faint supernovae. Therefore, we can emanine the impact of faint supernovae on the abundance of certain elements. We have made the approximation that the yields of 20MSun pre-supernovae stars are representative of faint supernovae in the 15-40 MSun mass range.

STELLAB (Stellar abundances) plots real observational data for different galaxies, and is plotted with OMEGA for comparison. SYGMA (Stellar Yields for Galactic Modelling Applications) calculates the chemical composition of the ejecta of Simple Stellar Populations (SSP) as a function of metallicity and time (Côté et al. 2016b) Chem_evol couples the stellar yields with the the distribution of initial stellar masses (Initial Mass Function [IMF]).


In [ ]:
#import modules
#sygma and omega share the same chem_evol class

import chem_evol
import sygma
import omega
import stellab

#import Python plotting packages
import matplotlib
import matplotlib.pyplot as plt

#Define Stellab
stellab = stellab.stellab()
%matplotlib nbagg

The Initial Mass Function

The Initial Mass Function describes the initial mass distribution of stars in a stellar population. Applying limits to this function allows us to control the mass range of the data it couples with. In this case it limits the range of the extra yields data. We have chosen the mass range 15-40 (Solar Mass), the mass range of the 1997D and 1999br faint supernovae Nomoto et al. 2013. Every star within this mass range is assumed to release the same ejecta as the 20MSun model. In reality, stars with different masses will create ejecta of different masses and composition. The 20 Msun star model is however an average representative of the 15-40MSun mass range. Z (iniZ) refers to the metallicity of the stellar population (one starburst).


In [ ]:
# Run a SYGMA instance to access the Initial Mass Function (IMF)
s_imf = sygma.sygma(iniZ=0.02, mgal=1.0)

# Define the mass range in which the extra yields will be applied
m_low = 15.0
m_up = 40.0

In [ ]:
# Calculate the number of stars in that mass range per units of Msun formed
A = IMF = 1.0 / s_imf._imf(s_imf.imf_bdys[0], s_imf.imf_bdys[1], 2)
nb_extra_star_per_m = A * s_imf._imf(m_low, m_up, 1)
print (nb_extra_star_per_m)

The model assumes that all the faint SNs occur between 7.581 x 10^6 and 9.588 x 10^6 years. These are the respective lifetimes of the 20MSun and 25MSun star models at Z = 0.02. The number of faint SN per unit solar mass is constant throughout this time period (~0.00368 / MSun). The t and R arrays define the extra source time-distribution rate for a simple stellar population. DTD stands for delay-time distribution. The extra pre-supernovae data yields are contained within 'yield_tables/extra_source_20Msun_O.txt'. The yield tables for preSN of different mass ranges are needed for future simulations for more precise analysis of faint SN ejecta.


In [ ]:
# Create the DTD and yields information for the extra source
# ==========================================================
# Event rate [yr^-1] as a function of time [yr].
# This assumes that all extra yields will be ejected
# between 7.581E+06 and 9.588E+06 years (the lifetimes
# of the 20 and 25 Msun models at Z = 0.02).
t = [7.580E+06, 7.581E+06, 9.588E+06, 9.589E+06]
R = [0.0, 1.0, 1.0, 0.0]

In [ ]:
# Build the input DTD array
dtd = []
for i in range(0,len(t)):
    dtd.append([t[i], R[i]])
    
# Add the DTD array in the delayed_extra_dtd array.
delayed_extra_dtd = [[dtd]]

In [ ]:
# Define the total number of event per unit of Msun formed.  
delayed_extra_dtd_norm = [[nb_extra_star_per_m]]

# Define the total mass ejected by an extra source
# Here, it would be best to find a correction factor 
# to account for the different total mass ejected by
# stars having different masses. For now, each star
# in the mass range eject the same ejecta as the 20Msun
# model.
delayed_extra_yields_norm = [[1.0]]

# Define the yields path for the extra source
extra_yields = ['yield_tables/extra_source_20Msun_O.txt']
delayed_extra_yields = extra_yields

Milky Way Commands

Omega incorporates the star formation history and data of different galaxies. The Milky Way is our chosen galaxy, although this can be changed using the galaxy parameter. The initial Milky Way mass is assumed to be 1e11 solar masses. The transition mass refers to the transition between AGB (Asymptopic Giant Branch)stars and massive stars. A transition mass of 8Msun indicates that all stars below 8Msun in the IMF will have AGB star yields attached to them, while the stars above 8 Msun will have massive star yields attached to them.


In [ ]:
# Use a different iPython notebook cell.
# OMEGA simulation with pre-SN yields and transition mass = 8 Msun
o = omega.omega(galaxy='milky_way', mgal=1e11, delayed_extra_dtd=delayed_extra_dtd, \
                delayed_extra_dtd_norm=delayed_extra_dtd_norm, delayed_extra_yields=delayed_extra_yields, \
                delayed_extra_yields_norm=delayed_extra_yields_norm, transitionmass=8.0)

# OMEGA simulation with pre-SN yields and transition mass = 10 Msun
o_trans = omega.omega(galaxy='milky_way', mgal=1e11, delayed_extra_dtd=delayed_extra_dtd, \
                delayed_extra_dtd_norm=delayed_extra_dtd_norm, delayed_extra_yields=delayed_extra_yields, \
                delayed_extra_yields_norm=delayed_extra_yields_norm, transitionmass=10.0)

# OMEGA simulation without pre-SN yields and transition mass = 8 Msun
o_no = omega.omega(galaxy='milky_way', mgal=1e11, transitionmass=8.0)

# OMEGA simulation without pre-SN yields and transition mass = 10 Msun
o_no_trans = omega.omega(galaxy='milky_way', mgal=1e11, transitionmass=10.0)

Plotting the stellar abundance ratios

We can now plot the stellar abundance ratios. The element(s) of interest can be changed by altering the xaxis/yaxis identities below. In this case, O/Fe is plotted against Fe/H. The abundance ratios are logarithmic as shown below. $$[A/B]=\log(n_A/n_B)-\log(n_A/n_B)_\odot.$$ where A and B are the relevant elements. This equation links $n_A/n_B$ (the abundance ratio for a star) and $(n_A/n_B)_\odot$ (the solar normalisation). This ensures the resulting ratios are notated relative to the Suns own abundance ratios. If [A/B]= 0 this indicates that the star shares the same element abundance ratio for A and B. A star therefore has an abundance ratio $10^{A/B}$ times that of the Sun.The solar normalisation can be extracted from a reference paper, using the norm parameter, which in this case is 'Grevesse_Noels_1993'.


In [ ]:
#Obtain a list of normalisation reference papers
stellab.list_solar_norm()

The observational data included by STELLAB can be changed by altering the list of references in obs. It is recommended that this list is kept relatively compact. Adding too much data may dominate the graph, making it difficult to observe important relationships. We have removed APOGEE, as it clouds particular areas of the graph, and does not currently feature notable error data. APOGEE can of course be returned if necessary. A list of references can be obtained using the command below.


In [ ]:
stellab.list_ref_papers()

The abundance ratios that are plotted can be changinging by swapping the element symbols in the x/y axis brackets. A list of included data is incorporated in obs, which STELLAB will plot. References can be removed from this list to reduce data dispersion and clouding. Error bars are included for detailed and accurate analysis of the ratios. These can be turned on by including 'show_err=True, show_mean_err=True' within the plot_spectro function above. However the inclusion of error may lead to OMEGA being plotted underneath STELLAB,so it is important to use an overplotting procedure as shown above. OMEGA plots in white underneath the coloured lines in this instance. This makes it easier to differentiate between OMEGA plots and STELLAB data. The colours are changed by changing the 'color' parameter. See Dale et.al for colour examples. Not all STELLAB data currently contain error bars, as many of these papers do not include them in yield tables, and they instead state an average error. These error bars can be plotted separately on the graph. The error markers respond to the STELLAB data which match their shape and colour.


In [ ]:
#Define your X and Y axis
xaxis = '[Fe/H]' 
yaxis = '[O/Fe]'
# Plot observational data
obs = [
'stellab_data/milky_way_data/Venn_et_al_2004_stellab',
'stellab_data/milky_way_data/Akerman_et_al_2004_stellab',
'stellab_data/milky_way_data/Andrievsky_et_al_2007_stellab',
'stellab_data/milky_way_data/Andrievsky_et_al_2008_stellab',
'stellab_data/milky_way_data/Andrievsky_et_al_2010_stellab',
'stellab_data/milky_way_data/Bensby_et_al_2005_stellab',
'stellab_data/milky_way_data/Bihain_et_al_2004_stellab',
'stellab_data/milky_way_data/Bonifacio_et_al_2009_stellab',
'stellab_data/milky_way_data/Caffau_et_al_2005_stellab',
'stellab_data/milky_way_data/Gratton_et_al_2003_stellab',
'stellab_data/milky_way_data/Israelian_et_al_2004_stellab',
'stellab_data/milky_way_data/Nissen_et_al_2007_stellab',
'stellab_data/milky_way_data/Reddy_et_al_2003_stellab',
'stellab_data/milky_way_data/Spite_et_al_2005_stellab',
'stellab_data/milky_way_data/Battistini_Bensby_2016_stellab',
'stellab_data/milky_way_data/Nissen_et_al_2014_stellab',
'stellab_data/milky_way_data/Ramirez_et_al_2013_stellab',
'stellab_data/milky_way_data/Bensby_et_al_2014_stellab',
'stellab_data/milky_way_data/Battistini_Bensby_2015_stellab',
'stellab_data/milky_way_data/Adibekyan_et_al_2012_stellab',
'stellab_data/milky_way_data/Aoki_Honda_2008_stellab',
'stellab_data/milky_way_data/Hansen_et_al_2012_pecu_excluded_stellab',
'stellab_data/milky_way_data/Ishigaki_et_al_2012_2013_stellab',
'stellab_data/milky_way_data/Roederer_et_al_2009_stellab',
'stellab_data/milky_way_data/Roederer_et_al_2014_pecu_excluded_stellab']



stellab.plot_spectro(xaxis=xaxis, yaxis=yaxis, galaxy= 'milky way', norm='Grevesse_Noels_1993', show_err=True, show_mean_err=True, obs=obs, ms=4)

# Extract numerical predictions
xy = o.plot_spectro(xaxis=xaxis, yaxis=yaxis, return_x_y=True)
xy_trans = o_trans.plot_spectro(xaxis=xaxis, yaxis=yaxis, return_x_y=True)
xy_no = o_no.plot_spectro(xaxis=xaxis, yaxis=yaxis, return_x_y=True)
xy_no_trans = o_no_trans.plot_spectro(xaxis=xaxis, yaxis=yaxis, return_x_y=True)

# Plot white lines - these make it easier to differentiate OMEGA plots from STELLAB data
plt.plot(xy[0], xy[1], color='w', linewidth=3, zorder=999)
plt.plot(xy_no[0], xy_no[1], color='w', linewidth=3,zorder=999)
plt.plot(xy_trans[0], xy_trans[1], color='w', linewidth=3, alpha=0.9, zorder= 999)
plt.plot(xy_no_trans[0], xy_no_trans[1], color='w', linewidth=3, alpha=0.9, zorder=999)

# Plot coloured lines
plt.plot(xy[0], xy[1], color='g', linewidth=1.5, label='PreSN for M=[15-40], M_trans=8',zorder=1000)
plt.plot(xy_trans[0], xy_trans[1], color='r', linewidth=1.5, linestyle='--', label='PreSN for M=[15-40], M_trans=10', zorder=1000)
plt.plot(xy_no[0], xy_no[1], color='midnightblue', linewidth=1.5, label='Original, M_trans=8', zorder=1000)
plt.plot(xy_no_trans[0], xy_no_trans[1], color='yellow', linewidth=1.5, linestyle='--', label='Original, M_trans=10', zorder=1000)

# Update the legend
plt.legend(loc='center left', bbox_to_anchor=(1.01, 0.5), markerscale=1, fontsize=10)

# x and y positions of the added error bars
xlow = -4.
xhigh = 0.5
ylow = -1.4
yhigh = 1.8
plt.xlim(xlow,xhigh)
plt.ylim(ylow,yhigh) 

#The size,shape and colour of the error bars and their markers
yerror=[0.16,0.05,0.0,0.23,0.2, 0.01] #dex
xerror=[0.07,0.03,0.03,0.05,0.24, 0.2] #size of error
marker=['bs','cx','bo','cs', 'co', 'rs'] #marker shape
markers=[4,4,4,4,4, 4] #marker size
colorbar=['b','c','b','c', 'c', 'r']
shift_x=[0.1,0.25, 0.37, 0.48, 0.8, 1.3] # these shift the position of the error bars
j=0
for i in yerror:
    plt.errorbar(xhigh-shift_x[j],yhigh-0.26,xerr=xerror[j],yerr=i,fmt=marker[j],markersize=markers[j],ecolor=colorbar[j],capsize=2)
    j=j+1

The impacts of the faint supernovae on the Milky Way can be observed by comparing the PreSN plots to the observational STELLAB plots and the Original (no PreSN) plots.


In [ ]: